rm(list = ls())
require(stats4)
require(nlWaldTest)
require("multiwayvcov")
require(msm)
require(stargazer)
require(mfx)

data21=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december6noon_type7clean.csv", header=TRUE)
data22=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_december15_type7clean.csv", header=TRUE)
data23=read.csv("C:/Users/naten/OneDrive - University of Tennessee/Dean Stuff/data/columbia_experiment_january26_type7(clean).csv", header=TRUE)


dataframe21=data.frame(data21)
dataframe22=data.frame(data22)
dataframe23=data.frame(data23)

#some definitions
expframe=rbind(dataframe21,dataframe22,dataframe23)
expframe=subset(expframe,uid>800)
expframe$correct <-(expframe$state==4&(expframe$response==2|expframe$response==7|expframe$response==8|expframe$response==9))|(expframe$state==3&(expframe$response==1|expframe$response==6|expframe$response==3|expframe$response==5))
expframe$incentive <- (expframe$qid==5)*5+(expframe$qid==1)*40+(expframe$qid==6)*70+(expframe$qid==7)*95
accuracy=c(sum((expframe$incentive==5)*expframe$correct)/sum(expframe$incentive==5),sum((expframe$incentive==40)*expframe$correct)/sum(expframe$incentive==40),sum((expframe$incentive==70)*expframe$correct)/sum(expframe$incentive==70),sum((expframe$incentive==95)*expframe$correct)/sum(expframe$incentive==95))
incentivec=c(5,40,70,95)
uidlist=unique(expframe$uid)

#check N
#**
length(uidlist)
#**

#for ease of factor use
expframe$incentive5=(expframe$incentive==5)
expframe$incentive40=(expframe$incentive==40)
expframe$incentive70=(expframe$incentive==70)
expframe$incentive95=(expframe$incentive==95)

#Check NIAC 
expframe$Bresponse <-(expframe$response==2|expframe$response==7|expframe$response==8|expframe$response==9)
expframe$Bstate <- expframe$state==4
expframe$Aresponse <-(expframe$response==1|expframe$response==6|expframe$response==3|expframe$response==5)
expframe$Astate <- expframe$state==3

#NIAC in aggregate

#Base Regression
expframe$incentivedec=expframe$incentive/100
modelNIAC=lm(Aresponse~factor(Astate+incentivedec)+0,data=expframe)

#Adjust COV
vcovmodNIAC=vcovHC(modelNIAC, type="HC3")
vcovmodNIACcluster=cluster.vcov(modelNIAC, expframe$uid)


Textstr=vector("character")
#Aggregate test
Textstr[1]="((b[6]+(1-b[2]))-(b[5]+(1-b[1])))"
waldNIAC1=nlWaldtest(modelNIAC,texts=Textstr,Vcov=vcovmodNIACcluster)
waldNIAC1

Textstr[1]="((b[7]+(1-b[3]))-(b[5]+(1-b[1])))"
waldNIAC2=nlWaldtest(modelNIAC,texts=Textstr,Vcov=vcovmodNIACcluster)
waldNIAC2

Textstr[1]="((b[8]+(1-b[4]))-(b[5]+(1-b[1])))"
waldNIAC3=nlWaldtest(modelNIAC,texts=Textstr,Vcov=vcovmodNIACcluster)
waldNIAC3

Textstr[1]="((b[7]+(1-b[3]))-(b[6]+(1-b[2])))"
waldNIAC4=nlWaldtest(modelNIAC,texts=Textstr,Vcov=vcovmodNIACcluster)
waldNIAC4

Textstr[1]="((b[8]+(1-b[4]))-(b[6]+(1-b[2])))"
waldNIAC5=nlWaldtest(modelNIAC,texts=Textstr,Vcov=vcovmodNIACcluster)
waldNIAC5

Textstr[1]="((b[8]+(1-b[4]))-(b[7]+(1-b[3])))"
waldNIAC6=nlWaldtest(modelNIAC,texts=Textstr,Vcov=vcovmodNIACcluster)
waldNIAC6


NIACOLSAggvec=c(waldNIAC1$p.value,waldNIAC2$p.value,waldNIAC3$p.value,waldNIAC4$p.value,waldNIAC5$p.value,waldNIAC6$p.value)

#TableA2.2 column 1******************
NIACOLSAggvec
#******************




#NIAC in aggregate for each pair to check for significant differences
NIACaggframe=data.frame(incentive_level_1=vector('numeric'),incentive_level_2=vector('numeric'),high_inc_avgmeffect=vector('numeric'),pvalue=vector('numeric'))
for(i in 1:3){
for(j in (i+1):4){
playframe=subset(expframe,incentive==incentivec[i]|incentive==incentivec[j])
playframe$highinc=playframe$incentive==incentivec[j]
playmodel=logitmfx(correct~highinc,data=playframe,robust=TRUE, clustervar1 = "uid",atmean=FALSE) 
stargazer(playmodel$mfxest)
stargazer(playmodel$fit)
addframe=data.frame(incentive_level_1=incentivec[i],incentive_level_2=incentivec[j],high_inc_avgmeffect= playmodel$mfxest[1,1], pvalue=playmodel$mfxest[1,4])

NIACaggframe=rbind(NIACaggframe,addframe)
print(i)
print(j)

}
}

#NIAC in aggregate for each pair to check for significant differences
NIACaggframe=data.frame(incentive_level_1=vector('numeric'),incentive_level_2=vector('numeric'),high_inc_avgmeffect=vector('numeric'),pvalue=vector('numeric'))
for(i in 1:3){
for(j in (i+1):4){
playframe=subset(expframe,incentive==incentivec[i]|incentive==incentivec[j])
playframe$highinc=playframe$incentive==incentivec[j]
playmodel=logitmfx(correct~highinc,data=playframe,robust=TRUE, clustervar1 = "uid",atmean=FALSE) 
stargazer(playmodel$mfxest)
stargazer(playmodel$fit)
addframe=data.frame(incentive_level_1=incentivec[i],incentive_level_2=incentivec[j],high_inc_avgmeffect= playmodel$mfxest[1,1], pvalue=playmodel$mfxest[1,4])

NIACaggframe=rbind(NIACaggframe,addframe)
print(i)
print(j)

}
}

#TableA2.2 column 2******************
stargazer(NIACaggframe,rownames=FALSE,summary=FALSE)
#******************
